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Abstract 

Post-transductional modifications tune the functions of proteins and regulate the 
collective dynamics of biochemical networks that determine how cells respond to 
environmental signals. For example, protein phosphorylation and nitrosylation are 
well-known to play a pivotal role in the intracellular transduction of activation and 
death signals. A protein can have multiple sites where chemical groups can reversibly 
attach in processes such as phosphorylation or nitrosylation. A microscopic descrip- 
tion of these processes must take into account the intrinsic probabilistic nature of 
the underlying reactions. We apply combinatorial considerations to standard en- 
zyme kinetics and in this way we extend to the dynamic regime a simplified version 
of the traditional models on the allosteric regulation of protein functions. We link 
a generic modification chain to a downstream Michaelis-Menten enzymatic reaction 
and we demonstrate numerically that this accounts both for thresholds and long 
time delays in the conversion of the substrate by the enzyme. The proposed mech- 
anism is stable and robust and the higher the number of modification sites, the 
greater the stability. We show that a high number of modification sites converts a 
fast reaction into a slow process, and the slowing down depends on the number of 
sites and may span many orders of magnitude; in this way multisite modification of 
proteins stands out as a general mechanism that allows the transfer of information 
from the very short time scales of enzyme reactions (milliseconds) to the long time 
scale of cell response (hours). 
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1 Introduction 



With the advancement of biochemical techniques a huge amount of data on the 
post-trasductional modification of proteins and on its role in the regulation 
of signal propagation through biochemical networks is now available. This 
knowledge is challenging our understanding of the cells' behaviour, and efforts 
by experts from different disciplines are now required to sort the basic dynamic 
principles out of experimental observations. As it has been recently pointed 
out [1], this process is not obvious at all because of cultural differences that 
drive scientists from different disciplines either to prefer the deep molecular 
details of biochemical networks or to tackle the generic underlying principles. 
Here we follow this latter perspective, and investigate the general dynamic 
consequences of multiple chemical modifications of proteins within biochemical 
networks 

Proteins can be modified by the attachment and detachment of various chemi- 
cal groups and by means of different mechanisms. Reversible chemical modifi- 
cations of proteins are now recognized to play a pivotal role in the regulation of 
the dynamic behaviour of biochemical networks, and therefore in the response 
of cells to environmental signals. The phosphorylation/dephosphorylation of 
tyrosine residues, for example, allows the propagation of signals from the cell 
surface to the nucleus thus enabling the transcription of specific genes in re- 
sponse to the presence of environmental activation molecules [2,3,4]; another 
example is the nitrosylation/denitrosylation of cysteine residues as a conse- 
quence of the redox state of the intracellular environment which can activate 
enzymes, such as caspase-3, that partecipate in the activation of the apoptotic 
program that leads to controlled cell death [5]. Cell activation, proliferation 
and death, in turn, are at the basis of animal physiology and pathology. For 
example, the control of the cell cycle through protein phosphorylation is im- 
portant for the activation of an immune response against foreign antigens 
[6]; the activation of the apoptotic program regulates tissue homeostasis and 
prevents the onset of autoimmune diseases and cancer [7]. 

Central to most biological networks is the existence of biochemical paths that 
behave dynamically as on-off switches [8,9]. The chemical modification of pro- 
teins on multiple aminoacid residues, like e.g. multisite phosphorylation, has 
been argued to confer a switch-like character to these proteins. One outstand- 
ing example is the retinoblastoma protein (Rb) with its 16 putative phospho- 
rylation sites, where at least 10 of them must be phosphorylated by cyclin- 
dependent kinases to promote the abrupt and irreversible Gl-S transition 
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along the cell cycle [10,11]. The switch-like behavior of proteins within bio- 
chemical networks has been traditionally modeled using the phenomenological 
Hill function [12]: 

k-X n 

= TTkX" (1) 

where X is the concentration of some biochemical species, A; is a positive 
constant and n is the so-called Hill coefficient. The Hill function provides 
switch-like sigmoidal behaviors for large n values and models phenomenolog- 
ically the multiple interaction of X with n other molecules. And yet, while 
the Hill function is quite good for phenomenological work it is well-known to 
be physically untenable for the following reasons: 1) the Hill equation implies 
simultaneous molecular interactions and this does not reflect a possible reac- 
tion scheme [13]; 2) in the real world, the concentration of proteins, enzymes 
and substrate fluctuates randomly as a consequence of molecular diffusion 
processes. In addition, molecules are randomly distributed between daughter 
cells at mitosis and this introduces an additional stochasticity in the molec- 
ular concentrations through cellular generations [14]: the Hill function does 
not incorporate the stochastic nature of the microscopic physical world and 
cannot be modified in this sense; 3) the sigmoidal Hill function becomes a step 
function, and thus describes a true threshold, only for n^oo, but experimen- 
tally, the chemical modification of just one key protein residue can turn the 
function of that protein on or off. For example, each subunit of the methion- 
ine adenosyl transferase (MAT) - a member of the caspase family of cysteine 
proteases - has 10 free cysteines but only cysteine 121 inhibits the activity of 
the enzyme when targeted by nitrosylation [5]. The Hill function with n — 1 
reduces to a standard Michaelis-Menten function which is no longer sigmoidal 
and therefore does not model a threshold behavior at all. The conclusion is 
that the Hill phenomenology must be replaced by deeper microscopic models 
to obtain a better understanding of the biochemical thresholds. 

From a chemical perspective, the processes leading to the chemical modifica- 
tion of proteins are equivalent to the allosteric regulation of enzymes whereby 
the substrate itself (homotropic interaction) or molecules other than the sub- 
strate (heterotropic interactions) can tune the activity of an enzyme. The 
allosteric theory dates back to the '60s when detailed models unifying both 
type of interactions where also developed, such as the famous models by 
Monod, Wyman and Changeaux (MWC) and by Koshland, Nemethy and 
Filmer (KNF) [15,16]. It is not the aim of the present work to discuss these 
models, their differences and their adherence to experimental data, as these 
aspects have been fully reviewed elsewhere [17,18]. Some general points should 
nonetheless be analyzed. Both the MWC and the KNF models deal with the 
microscopic interactions between a chemical species and the subunits of an 
enzyme, the chemical species being capable of stabilizing (MWC) or induc- 
ing (KNF) a conformational change in the quaternary structure of the target 
subunit (s). The conformational change of the quaternary structure - but the 
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theory has been recently extended to the conformational change of the ter- 
tiary structure of an enzyme within the very same framework of the MWC 
model [19] - influences in turn the enzymatic activity. In other words, these 
models focus on the mechanics (i.e. modifications of the protein structure) 
and not on the dynamics of the allosteric effect. And in fact, the models rely 
on a set of assumptions whose common one is the steady state hypothesis: 
protein complexes are in equilibrium with a very large number of ligands 
whose concentration is always constant [15,16,17,18]. This assumption allows 
the authors to linearize the model equations which, however, do still contain 
a number of parameters whose values must be estimated by fitting equations 
to experimental data. It has been pointed out that the high number of ad- 
justable parameters hampers the validation of a given allosteric model [17] 
and, we add, also hinders the inclusion of such equations into large biochemi- 
cal ensembles. Simpler models that do take into account the real variations of 
the ligand concentrations within a cell are therefore required to describe the 
dynamic behaviour of protein networks. Of course a model of this kind must 
be based on a simplification of the molecular details underlying the allosteric 
effect. We propose here to focus the analysis on the dynamic properties of 
the allosteric effect and we put forward a simple and very general microscopic 
model of the multiple chemical modifications of proteins. Our analysis does 
not explain how the allosteric effect is realized structurally by a given protein 
but it shows how certain important properties of biochemical networks - such 
as thresholds, long delays and stability - emerge from a simplified analysis of 
allosteric dynamics. 



2 Equilibrium concentrations of the modified protein forms 

We start by considering the reaction scheme shown in figure 1, which is ba- 
sically equivalent to that considered by Monod et al. [15], where a chemical 
species B can modify a number of sites on protein A. We also assume that 
the sites are all equivalent and that the modification dynamics for each site 
is independent from those of the other sites: this means that we consider the 
states A n and the transition chain shown in figure 1. If the single chemical 
modification dynamics is fast with respect to the observation time we can 
forget the dynamics of the transition chain and concentrate instead on the 
equilibrium probabilities. Indeed, if p is the probability that any given site is 
chemically modified, the equilibrium probability P n that at least n molecules 
of the species B bind to A is given by the sum of the corresponding proba- 
bilities from a binomial distribution and therefore if the volume of reaction 
V contains va = NaV [A] molecules of type A, where Na is the Avogadro 
constant and the square brackets denote molar concentrations, the average 
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number of A molecules with at least n occupied sites is: 



N fN\ 

l=n \ I 



N-l 



(2) 



Every single modification reaction can be represented as a generic bimolecular 
reaction, A + B ^ AB, so that the following equations hold 



d[A] 

dt 
d[B]_ 

dt 
d [AB] 

dt 



-k+ [A] [B] + k_ [AB] 
-k+ [A] [B] + A:_ [AB] 
= k + [A] [B] - k_ [AB] 



(3) 



The reaction dynamics (3) leads to the conservation equations [A] + [AB] = 
[A] and [B] + [AB] = [B] , where [A] and [B] are the initial concentra- 
tions, while at equilibrium the derivatives vanish and the condition [AB] = 
(k + /k-)[A][B] holds; combining these conditions we find an equation that can 
easily be solved, and we obtain 



[AB] = U([A] + [B] + ^ 



A 



([A]o-[B]o) 2 + 2([A] +[B] )^-+(^ 



(4) 



Referring to figure 1 we reinterpret the result (4) as a concentration of occupied 
sites, so that the probability p is 



[occupied sites] 
N[Ajo 



A 



(N[A] - [B] y + 2(N[A] + [B] ) k Y + 



(5) 



(6) 



where we note in passing that the dissociation constant k_/k + depends on 
temperature (in general k_/k + = exp(AG/RT), and AG is the Gibbs free 
energy. 
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3 Threshold and threshold robustness 



Now we assume that the molecule A somehow switches to an activated form 
A* when at least n modified sites are occupied (see figure 1), and that the 
process is reversible. This might occur because of a conformational switch [20] 
or because A releases some reaction product, but the detailed description of 
this step is not important, what really matters is that a change occurs only 
if at least n sites are involved. It is worth noting that this analysis includes 
the case n — 1. Only few molecules are actually activated by this mechanism 
when p is small. For example, in a spherical cell with a radius of 5 /im and a 
corresponding volume ps 5- 10~ 16 m 3 , where proteins have concentrations below 
10 /xM, there are fewer than 3 • 10 3 molecules of each kind of protein, and even 
with the not so small value p ~ 0.1 we find, e.g., P n ~ 4.5 • 10~ 7 for N = 16 
and n — 10 (as for the Rb protein), so that in this case there is on average 
much less than 1 activated molecule. This means that molecular discreteness 
is important, and sets a natural threshold for the process set in motion by the 
activated form of A: the process cannot proceed at all if, on average, there is 
less than 1 activated molecule and the threshold is given by the probability 
p (and by the corresponding concentrations in equation (5)) that solves the 
equation v^Pn = 1- Since p < 1, equation (2) can be approximated by the 
simple power law 

v A P n *N A V[A] ^Jp n (7) 

and this means that for n = 10 as in the example, the threshold is very sharply 
defined (see figure 3). The threshold defined by the solution of the equation 
v aPu — 1 j with n > 1, is not only sharp but very robust as well. There are 
many sources of potentially destructive fluctuations, like the variable number 
of A molecules that could change as a consequence of the unequal distribution 
of enzymes and substrates between daughter cells at mitosis [14] , but it is easy 
to see that these fluctuations are damped down, and that the damping-off is 
more successful for large n's. When the k_/k + ratio is negligible with respect 
to [A] and [B] , the probability (5) can be approximated as follows: 



o 



[A]o 



[B] Q /(N[A] ) if N > [B] /[A] 

1 if N < [B] /[A] 
(8) 

with a breakpoint at [B] = N[A] . In general, for low concentration [B] , 
expression (5) depends linearly on [B] : 
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Using the approximate expression (9) and equation (7) we find the threshold 
value: 

[B]o,thr = (n[A} + ^ [n a V[A] Q | (10) 

If we use standard formulas for the analysis of random fluctuations of experi- 
mental quantities, we find 

var[ff] , tfer 1 vaiV ( N[A] (n - 1) - (fc_/fc+) \ 2 var[A] 
([S] 0>thl .) 2 n 2 ^ 2 { n(N[A] + (k_/k + )) ) (L4] ) 2 1 ' 

and the end result is that volume fluctuations are damped by a factor n, 
so that, e.g., a 10% volume fluctuation produces a mere 1% change in the 
threshold position for n — 10, while the damping of substrate fluctuations 
varies from slightly sublinear (so that a 10% fluctuation of [A} produces a 9% 
change in the threshold position for n = 10) to a strong damping condition 
like in the case of volume (for k^/k + 3> iV[A| ) . 



4 Downstream Michaelis-Menten step 



As shown in figure 1, molecules A activated by multiple chemical modifica- 
tions are supposed to release an enzyme E that catalyses a Michaelis-Menten 
reaction where a substrate S is converted into a product R. The downstream 
reaction catalyzed by E writes: 

E + sh ES E + R 

and we use the Michaelis-Menten equations with the usual quasi-steady state 
assumption (QSSA, see, e.g. ref. [21] and references therein), so that the sub- 
strate consumption is described by the single differential equation 

d[S] ^ k 3 h[E] [S] k 3 [E] Q [S] 

dt ~ ki [S] + (k 2 + h) K m +[S] { ] 

where K m = (k 2 +ks) / k\ and the production of R is related to the consumption 
rate by the equality d[R]/dt = —d[S]/dt From the considerations above, it is 
easy to see that the concentration of the enzyme E equals the concentration 
of the molecules A with at least n occupied sites: 

[E]o = [A] j:{ N i y i (l-p) N - 1 (13) 

then, if we combine equations (12) and (13), we can integrate numerically 
equation (12). Figure 4 shows the behaviour of [R] obtained with parameters 
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in a range common to many biochemical reactions in the cell. Here we took 
the concentration [B] to be a linear time-dependent quantity: [B] = 16 ■ 
(10/iM) • (10 _5 s -1 )t, so that [B] reaches the critical value N[A] at t = 10 5 
s. Figure 5 shows that as the chemical modification level n grows the time 
delay At (time taken to reach 50% of the maximum product concentration) 
spans several orders of magnitude (from a few seconds up to a few hours). The 
production of R is not just delayed, it is also slowed down: figure 5 shows that 
the production interval spans several orders of magnitude as well. 



5 Conclusion 

In the mathematical modeling of cooperative reaction kinetics, threshold ef- 
fects are commonly achieved with phenomenological Hill equations, however 
here we no longer require this kind of phenomenological modeling and both 
thresholds and time delays have a clear meaning. This is particularly im- 
portant in models of the behavior of complex biochemical systems, as time 
delays coupled to negative feedback regulatory circuits are required to pro- 
duce bistability and limit cycle behavior [22]. We have neglected the coupling 
between the modification step and the downhill Michaelis-Menten step, be- 
cause the probabilistic model does not describe an actual reaction dynamics. 
The Michaelis-Menten step assumes that part of the enzyme is bound with the 
substrate to form the complex ES, and this acts as a store of E, which is re- 
leased with a characteristic relaxation time r = [ki[S] + (k 2 + ^3)] _1 [18]. For 
this reason we expect the downhill reaction to damp off any noise in the [B] 
concentration signal: thus the combination of the modification step with the 
downhill (irreversible) Michaelis-Menten step, makes the whole system even 
more robust, and resistant to environmental changes. 

All this is particularly intriguing from an evolutionary point of view: proteins 
with a high number of residues that can be modified chemically might in 
fact have been selected because this would intrinsically stabilize the threshold 
effect and hence the overall dynamics of the biochemical network operating in 
highly perturbed environments. 
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Fig. 1. a. Molecule A has N phosphorylation sites: here the process is schematically 
represented by molecules B that react with each site with on-off rates k+, b. 
When at least n sites out of the possible N sites are phosphorylated, molecule A 
releases an enzyme E which catalyses a Michaelis-Menten reaction that converts a 
substrate S into a product R. 




Fig. 2. Schematic representation of the phosphorylation-dephosphorylation chain. 
We assume that all sites are equivalent and that the enzyme state depends only on 
the number of phosphorylated sites, so that there are only A + 1 different states 
which are labeled Ai. The arrows represent the phosphorylating-dephosphorylating 
transitions. 
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2 5 10 20 50 100 200 

Concentration [B] Q (piM) 

Fig. 3. Log-log plot of the (average) number of active A* molecules vs. [B], obtained 
from equations (2) and (5). The number of active A* molecules with at least n = 10 
phosphorylated sites over N = 16 possible sites has been calculated using equation 
(2) with the following parameters: [A]q = 10/uM; k-/k+ = 10~ 6 M (this is in the 
observed range for the dissociation constant of ATP with cyclin-dependent kinases 
[23]) and with a cell volume « 5 • 10~ 16 m 3 , so that the cell contains approximately 
3000 A molecules as in the example discussed in the text. The dashed line marks 
the threshold level vaP u = 1- 
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Fig. 4. a) Relative concentration [i?]/[5]o obtained from the numerical integra- 
tion of equation (12) with condition (13); in this case K m = ImM, k% = 10 4 s _1 , 
[A] = 10/xM, [S] = ImM, k^/k + = 10" 6 M, N = 16, n = 10, and with a linear 
growth law [B] = 16(10/zM) • (lO^s -1 )*: in this way the concentration of [B]q 
equals the critical concentration A/"L4] at t = 10 5 s. At is the time delay needed to 
reach 50% of the maximum product concentration [i?] maa; = [<S]o- b) Scaled deriva- 
tive (jw^j '■ tne width T (measured at 10% of the peak value) gives an estimate 
of the process duration. 
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Fig. 5. Plot of the time delay At and of the width T vs. the number of phospho- 
rylated sites n from the integration in figure 4. As n grows both At and V span 
several orders of magnitude, and range from a few seconds to several hours. 
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